Spatial tomography of light resolved in time, spectrum, and polarisation

Measuring polarisation, spectrum, temporal dynamics, and spatial complex amplitude of optical beams is essential to studying phenomena in laser dynamics, telecommunications and nonlinear optics. Current characterisation techniques apply in limited contexts. Non-interferometric methods struggle to distinguish spatial phase, while phase-sensitive approaches necessitate either an auxiliary reference source or a self-reference, neither of which is universally available. Deciphering complex wavefronts of multiple co-propagating incoherent fields remains particularly challenging. We harness principles of spatial state tomography to circumvent these limitations and measure a complete description of an unknown beam as a set of spectrally, temporally, and polarisation resolved spatial state density matrices. Each density matrix slice resolves the spatial complex amplitude of multiple mutually incoherent fields, which over several slices reveals the spectral or temporal evolution of these fields even when fields spectrally or temporally overlap. We demonstrate these features by characterising the spatiotemporal and spatiospectral output of a vertical-cavity surface-emitting laser.

I nformation can be encoded in optical beams by design, either in space through sculpted complex wavefronts [1][2][3] or, in time, through tailored pulse envelopes 4 , with further multiplexing of information channels possible using the spectral and polarisation dimensions. Encoded information can also be a consequence of natural phenomena such as the nonlinear interaction of light with the gain medium of the laser cavity [5][6][7] . In both scenarios, optical beams can be generally composed of multiple incoherent spatial fields, with variable time delays and with numerous spectral peaks. It is particularly challenging to fully characterise the spatial phase of the beam in this most general casethat is, where the beam is composed of multiple spatial components that may or may not be mutually coherent, may or may not exist at the same time, and may or may not have the same wavelength. This diversity of components carried by the beam makes decoding the complete information extremely difficult in the expanding frontier of applications ranging from imaging 8 to linear and nonlinear spatio-temporal beam-shaping [9][10][11] .
Conventional optical beam analysis techniques with the ability to recover not only the spatial intensity 12,13 but also spatial phase and coherence either utilise an external reference to recover the optical signal characteristics [14][15][16] or work on a self-referencing principle 17,18 . In the external reference case, the spatial amplitude and phase of light is measured by interfering the unknown optical beam with an external local oscillator, which imposes limitations on the spectral band and bandwidth that can be probed, with a suitable high-quality local oscillator not always available in the spectral region of interest. The self-referencing approaches, such as modal analysis 19,20 , alleviate the need for a local oscillator by using part of the beam itself as a reference. In this case, the spatial components of the unknown beam are determined via a series of projective intensity measurements using spatial correlation filters. The filters interfere probed modal components with a preselected internal reference mode to recover complex superposition coefficients of all spatial modes that constitute the original beam. The modal analysis method allows profiling of singlewavelength wavefronts 21,22 and dynamic tracking of temporal mode instabilities in lasers at a camera-speed [23][24][25] ; however, only if at least one mode is known a priori to exist in the beam 17 . Often, there is no such prior knowledge, or indeed the beam has no single mode that can serve as an appropriate reference. For example, even for omnipresent optical sources, such as laser diodes, the spatial components of the beam typically occupy spectrally discrete and distinct locations, which means there is no spatial component that can be singled out as a suitable reference over the whole spectral range. Similarly, an optical beam modulated in time has a dynamic wavefront composed of various spatial modes, with none guaranteed to exist over the whole temporal range and therefore none suitable as a reference. This severely narrows the range of applicability of self-referencing techniques and renders them incapable to analyse arbitrary optical beams.
A conceptually similar approach to modal analysis that can measure a completely arbitrary state, either mixed or pure, exists in quantum optics. In the context of spatial optical fields, a mixed state represents multiple mutually spatially incoherent fields whereas a pure state describes a single coherent field. By obtaining a sequence of projective tomographic measurements that are sensitive to certain physical aspects of the state, such as spatial features or polarisation, a complete representation of the state can be recovered in the form of the density matrix 26,27 . The density matrix formalism is applied to study a myriad of quantum states, ranging from vibrational modes in molecules 28 to an analysis of quantum gates 29 . A number of studies also used the technique to study the classical degrees of freedom of light individually, including in the spatial [30][31][32][33][34][35][36] , polarisation 37,38 and the time-frequency domain 39 and were even applied to hybrid entangled states with multiple degrees of freedom (space, polarisation) and (space, polarisation, time energy) with each degree of freedom probed by means of state tomography framework 40,41 . However, none of the existing methods provides access to polarisation resolved, high-dimensional spatial state density matrices of the unknown beam at its numerous spectral and temporal slices and therefore cannot reveal the rich, complex dynamics of arbitrary optical fields.
In this work, we use a high-dimensional (N = 21) spatial analogue of Stokes polarimetry 32,38 (Fig. 1a, b) to reconstruct an unknown optical beam in space along with its temporal, spectral and polarisation characteristics via a series of spatial projective measurements. Using polarisation optics, together with an oscilloscope and a spectrometer as illustrated in Fig. 2a, we typically acquire 2 × N t + 2 × N s spatial state density matrices in a highly parallel fashion for the two polarisations, N t = 1000 temporal and N s = 250 spectral slices. The spatial dimensions (N), as well as the number of temporal (N t ) and spectral (N s ) bins, can be readily increased if the application requires it. We perform the spatial projective measurements using a tandem of a spatial light modulator (SLM) and a single-pixel detector in the form of singlemode fibre (SMF) as illustrated in Fig. 1b. Each projection measurement corresponds to a pairwise interference between two potential spatial components (modes) of the unknown beam, with the full measurement interfering all the possible spatial components that exist within the beam.

Results
In stark contrast to the modal analysis approach that necessitates prior knowledge of a reference mode, our approach does not require a known reference state, and handles even scenarios with multiple references or no reference at all. Any spatial state that exists in any wavelength or temporal bin can be leveraged by Stokes method as a reference. If there is a mixture of incoherent states in a given bin, all the states are acting as independent references for all the other states. As a result, the Stokes method can spatially analyse the beam along the spectral and temporal axis, as long as there is detectable light in a given wavelength or temporal bin, and the spatial projective measurements span the existing spatial states in that bin.
The collection of light with the single-pixel detector in the form of SMF is modular in terms of routing and splitting the light to specialised detectors that are optimised for temporal and spectral resolution, such as oscilloscopes and spectrometers. This provides significant advantages over camera-based approaches that have limited refresh rates to resolve ultra-fast temporal dynamics, and low pixel counts to provide sufficient number of spectral bins. Probing the whole transversal profile of the beam at once by a single-pixel detector also improves the detection efficiency and decreases noise compared to raster-scanning or camera-based approaches that only examine a limited subspace of the beam at any time 42 . Interestingly, despite no spatial resolution of the detector, the spatial resolution of the measured field can be arbitrarily high, since the basis used for projections is numerically calculated.
We illustrate the merit of the technique by analysing the beam of a modulated vertical-cavity surface-emitting laser (VCSEL) diode which supports multiple mutually incoherent spatial modes within the beam, and can be readily modulated in time. The rich spatiotemporal and spectral output of the VCSEL resists full analysis using existing techniques despite its popularity in diverse products such as LiDAR for autonomous vehicles, face recognition and datacoms. Moreover, VCSEL attracts considerable attention in high-power applications, mainly from the perspective of fundamental studies such as spatio-temporal instabilities on a sub-nanosecond timescale 7 and wave-chaotic behaviour due to symmetry breaking in the laser cavity 5 .
Principle of projective measurements. In standard Stokes polarimetry illustrated in Fig. 1a, the unknown polarisation state of light, mathematically described by a density matrixρ in , is interrogated by a sequence of polarisation analyser states that are defined by the orientation of the polarisation optics. In the depicted example, the power coupled into the SMF is equal to the expectation value of the unknown input state (ρ in ) in the horizontal polarisation state. This horizontal polarisation state would exist at the input plane if the light was hypothetically launched from the SMF. Detecting the power with an SMF for multiple analyser states enables the recovery of the Stokes vector (S), with its elements S n signifying weights of Pauli matricesσ n in a sum that reconstructs the density matrix (ρ measured ) of the unknown state. The density matrix encapsulates complete information about the unknown input state, including whether it represents a mixture of up to N = 2 states or a pure state. The situation is mathematically equivalent for spatial states (Fig. 1b), however, to perform the projective measurements experimentally in this case, we have to substitute the polarisation optics for a spatial light modulator, which acts as a reconfigurable spatial analyser state device. The power coupled into the SMF in this scenario is equal to the expectation value (|field overlap| 2 ) of the unknown input (ρ in ) in the state defined by LP11 mode with lobes oriented horizontally, as illustrated in Fig. 1b. We can again imagine that we hypothetically launch the light from the SMF towards the SLM which transforms the fundamental LP01 mode of the SMF into LP11 projective state at the source plane. We can mathematically overlap the source and the hypothetical LP11 field at the source plane but since no detector exists there the overlap value cannot be determined experimentally at this plane. However, the reciprocity of light ensures conservation of field overlap value in the system, which in turn allows its experimental recovery at the SMF plane (more details in Supplementary Note 4).
Experimental setup. Figure 2 (a) shows a simplified schematic of the experimental setup. The unknown, pure or mixed spatial state of light that generally depends on wavelength and time (ρ in λ; t ð Þ) is separated by a polarisation beam splitter (PBS) into two paths to allow independent interrogation of vertical (V) and horizontal (H) polarisation by a single, large-area SLM that handles both polarisations. The orthogonal polarisation paths are recombined with another PBS and coupled into an SMF. The SMF subsequently routes the light to an optical spectrum analyser (OSA) and oscilloscope to analyse the input spatial state with respect to wavelength and temporal degrees of freedom. Identical temporally modulated copies of the beam are generated for each SLM projection mask via the pseudo-random binary sequence (PRBS) generator that modulates the laser in a synchronised manner with the oscilloscope. More detailed experimental setup configuration can be found in Supplementary Note 1, along with the calibration and alignment procedures of the system in Supplementary Note 2.
The unknown spatial state, represented by the density matrix (ρ in λ; t ð Þ), can have a spatial profile with fine features due to the existence of multiple coherent and incoherent spatial fields at each wavelength and time. In order to extend the technique to more complicated spatial fields with more features, necessitates the extension of Stokes formalism from the limited twodimensional case (Fig. 1b), where the analyser states are only sensitive to two modal components, into much higher dimension, Fig. 1 Working principle of Stokes polarimetry and its spatial analogue. An unknown pure or mixed state, described by density matrixρ in , passes through a sequence of polarisition (a) or spatial (b) analyser states encoded on a spatial light modulator (SLM). Analyser states are determined by eigenvectors of Pauli matricesσ n that represent all observables of a two dimensional Hilbert space. The intensity (expectation value) of a given observable is registered by a detector and a Stokes vector (S) is reconstructed by weighting the measured intensities with eigenvalues κ m n associated with eigenvectors ofσ n . The sum of Pauli matrices weighted by the Stokes vector elements (S n ) determines the measured density matrix. Measured mixed state can be graphically depicted within the volume of the Bloch sphere, whereas pure states reside on its surface 33 . until the Stokes states span all the modal components present in the beam. In our experimental case, a Stokes space dimension of N = 21 satisfies this condition, and actually exceeds it. When no a priori information about the beam is known, the Stokes dimensionality N should be chosen with some redundancy, at the price of carrying more projective measurements.
The natural generalisation of Pauli matrices, used for the reconstruction of the density matrix in two-dimensional scenario, into higher dimensions is facilitated by the generalised Gell-Mann matrices 43 that span the space of observables of an N-dimensional Hilbert space. A vector in such an N-dimensional Hilbert space, specifically its vector elements, can represent complex superposition coefficients of an optical field in any orthogonal spatial basis of interest. The choice of this basis can be arbitrary. To demonstrate that the tomographic measurement can be performed using a random orthogonal spatial basis, we generate the spatial analyser states in the following fashion. We take as a starting point a Laguerre-Gaussian (LG) basis (Fig. 2b), with the waist of w 0 = 3 μm reflecting the approximate extent of the investigated beam, and generate a single new spatial state composed of a random complex superposition of N states in the LG basis (Fig. 2c). We subsequently apply the Gram-Schmidt orthogonalisation process to create the remaining N − 1 orthogonal random spatial modes (Fig. 2c) to form the full random Laguerre-Gauss (RLG) basis with the dimension N.
The number of modes N in the initial selected basis should reflect the expected level of spatial details in the probed system. For instance, the spatial analyser states can sense high spatial frequency components only if the initial Stokes basis and the related random basis contains such high spatial frequency components. The choice of LG basis is convenient in cases with expected rotational symmetry, but the method is not limited to LG basis and works with Hermite-Gauss, position, Bessel, Hadamard or Fourier basis as per application needs. ð Þ is split into two optical paths by a PBS to allow independent interrogation of horizontal and vertical polarisation of the input state by the SLM. The recombined light is coupled into an SMF and guided into an oscilloscope and OSA for the analysis of the temporal and spectral degree of freedom. PBS -polarisation beam splitter, SLM -spatial light modulator, λ/2 -half-waveplate, SMF -single-mode fibre, OSA -optical spectrum analyser, PRBS -pseudo-random binary sequence generator. A complex superposition of Laguerre-Gauss (LG) spatial states (b), given by random complex coefficients c j , is used to construct one random Laguerre-Gauss (RLG) spatial state (c). The Gram-Schmidt orthogonalisation process generates the remaining states in the RLG basis (c). The k-th spatial analyser state (d) is constructed as a complex superposition of RLG spatial states with the complex superposition coefficients given by the k-th unique eigenvector of the Gell-Mann matrix set.
To generate the spatial analyser states that perform the tomographic measurements, we first find a set of unique eigenvectors of the Gell-Mann matrices. The vector elements of each eigenvector represent the complex superposition coefficients in the RLG basis and make up a given spatial analyser state (Fig. 2d). By sequentially applying the spatial analyser states using the SLM and measuring the intensity as a function of wavelength and time, using the OSA and oscilloscope respectively, we can reconstruct the generalised Stokes vector as both a function of wavelength and time. A blueprint for the calculation of 2N 2 − N = 861 analyser masks for Stokes dimension of N = 21 that generate the projective states, along with the mask alignment requirements is in Supplementary Note 2. Similar to the formalism presented in Fig. 1b, by summing the Gell-Mann matrices with weights given by the Stokes vector, we reconstruct the spectrally and temporally resolved density matrices of the unknown beam. The eigenvectors and eigenvalues of the reconstructed density matrices (expressed in the RLG basis), are the spatial eigenstates and their probabilities in the investigated wavelength or temporal bin. Since both the spectrum and the time trace are collected in a single sweep for one specific spatial analyser state, only 2N 2 − N measurements in total are required for the full reconstruction of the spatial, spectral and temporal information of the beam. A detailed theory of the highdimensional Stokes analysis is presented in Supplementary Note 3.
Spatio-spectral analysis. We chose the vertical-cavity surfaceemitting laser (VCSEL) diode as the unknown stateρ in λ; t ð Þ to demonstrate the method. VCSEL is a spatially, spectrally and temporally interesting source that represents a class of optical beams difficult to analyse using existing techniques. VCSEL cavity supports several spatial modes, with distributed spectral peaks (Fig. 3a). The position of the spectral peaks and the number of spatial modes varies based on the VCSEL driving conditions like the bias current and the temperature (see Supplementary Note 7 for modally resolved light-current curves and details of driving conditions). The number of spatial modes and their relative power impacts the overall spatial profile of the VCSEL, which influences its performance in applications such as face recognition and 3D scene scanning for augmented reality 44 . The ability to measure the spatio-spectral profile can provide useful feedback during VCSEL engineering and fabrication to improve the VCSEL performance in those applications.
In this section, we use the measured spectrally resolved spatial density matrix to distinguish the amplitude and phase of each spatial mode for each spectral slice. The spatial state density matrix unlocks the ability to resolve multiple mutually incoherent modes within the same spectral bin, which typically occurs in situations when spectral peaks of modes are closely packed and thus beyond the resolution of the spectrometer. We confirm the validity of the results by comparing the probability-weighted sum of reconstructed mode intensities with the intensity profiles obtained from a simple raster-scan, with results obtained very similar in both cases.
To perform the spatio-spectral analysis, we first collect the spectrum for each spatial analyser state in Fig. 2d and find the Stokes vector and the density matrix for each wavelength bin (resolution bandwidth of 70 pm given by the spectrometer) in the RLG basis. The measured density matrix at selected spectral peaks is depicted in Fig. 3b. The eigenvalues and eigenvectors of the density matrix represent the probabilities and the spatial states at each wavelength, respectively. We note that the acquired density matrix should theoretically be positive-semidefinite to ensure positive probability for each supported state. However, due to measurement-induced noise, this is not generally the case, and we have to apply Higham algorithm that finds the nearest positive semidefinite matrix in the Frobenius norm sense 45 .
The mutually incoherent spatial states with the three highest probabilities at each spectral peak of Fig. 3a are shown in Fig. 3c-e. The spatial states were constructed from the eigenvectors of the density matrix, with eigenvectors representing the complex superposition coefficients in the RLG basis. Despite the randomness of the RLG spatial basis, the reconstructed modes have the well-known linearly polarised (LP) mode spatial profiles typical for a circular laser cavity. Some spectral peaks have a significant second dominant state (Fig. 3d) and another weaker third state (Fig. 3e). This is most pronounced for the spatial modes corresponding to LP02 (green column), LP21 (red column) and LP11 (purple column) but also faintly visible for LP31 (orange column). In the LP11 case, the second dominant state is present simply due to leakage of the powerful neighbouring LP11 mode with a perpendicular orientation of lobes (brown). In the LP21 case, the wavelength bin contains a mixture of mutually incoherent states that are spectrally very close, within the 70 pm resolution of our spectrometer. Remarkably, even though we cannot spectrally resolve the modes, we can utilise the density matrix to recover not only the individual spatial profiles in the mixture but also their relative power. The ability of the technique to resolve spectrally overlapping spatial modes is one of its hallmark features. On the basis of mutual incoherence, it enables spatial differentiation of features beyond the spectral resolution. We explore this compelling attribute of Stokes analysis in detail in Supplementary Note 8, where we intentionally decrease the resolution bandwidth of the OSA to 10 nm, fully spectrally mixing the spatial modes in a single detected wavelength bin, yet the Stokes analysis yields the correct spatial states and their relative powers in the mixture. For the LP21 case, the recovered probabilities in the mixture are 63% for the LP21a and 27% for the LP21b state. Apart from the two dominant modes totalling 90% of power in the investigated wavelength bin, there is also approximately 2% of LP02 mode as visible on a partially distorted spatial profile in Fig. 3e. The 8% of unaccounted power distributes among the remaining 18 spatial eigenstates (not displayed) of the density matrix, all of them with a random, speckle-like spatial profile. There is no physical reason to expect that such random spatial states exist in the laser cavity. We instead associate the 8% of the remaining power to noise that inevitably occurs during the hundreds of sequential projections, either due to power fluctuation of the VCSEL or the thermally induced spectral offset.
The noise not only affects the amount of light that cannot be attributed to physical spatial modes but also influences the reconstruction fidelity of the spatial profile of the relatively weak modes. This is noticeable for the weak LP02 (Fig. 3e, red column) and also evident in the green column where both Fig. 3d, e are likely the weak, leaking LP21 modes. Nevertheless, despite these noise-induced effects, the mixed state reconstruction in Fig. 3f, obtained as a probability-weighted intensity sum of all 21 spatial eigenstates, is in almost perfect match with the intensity profiles (Fig. 3g) obtained by raster scan, which implies accurate determination of the spatial modes and their relative powers. We perform the raster scan by deflecting the beam using linear phase ramps (tilt) at the SLM, which scans the beam around the fixed core of the single-mode fibre. For each tilt, we collect the corresponding spectrum to generate the (x, y, λ) datacube, with relevant wavelength slices plotted in Fig. 3g. The influence of the noise level on the Stokes analysis is explored in-depth in Supplementary Note 9 for different levels of noise, RLG basis dimensionality, for pure and mixed states and also exploring the effect of orthogonality of input spatial states.  3 Spatio-spectral analysis of VCSEL (H polarisation). a The spectrum has multiple spectral peaks, each corresponding to one or more spatial modes. b The measured density matrix expressed in RLG basis for colour-coded spectral peaks. c-e The first three most dominant spatial eigenstates and their corresponding probabilities as reconstructed from the density matrix (b) using the RLG basis. f The probability-weighted intensity sum of all reconstructed eigenstates obtained via state tomography perfectly matches the intensity profiles (g) obtained by raster scanning the VCSEL beam over SMF by adding tilt on the SLM and collecting the spectrum for each tilt value. Quantitative agreement between (f) (tomography approach) and (g) (verification) is highlighted by the overlap integral value o evaluated in each spectral slice. The measured density matrix along with the reconstructed spatial eigenstates for each wavelength is in Supplementary Media 1 Spatio-temporal analysis. The number of spatial modes, their spectral peak position and relative power depends on the VCSEL temperature and bias. The applied bias does not have to be constant. VCSELs are often modulated at speeds of tens of GHz in short-reach interconnects that form the backbone of modern data centers and local area networks. The quality and speed of modulation can be influenced by mode competition within the cavity and the dynamic modal content of the beam impacts the data transmission through multimode fibres due to modal and chromatic dispersion. Having access to the temporal dynamics of VCSELs on an individual spatial mode basis is therefore of paramount importance in the quest to design faster optical interconnects.
The measured temporally resolved density matrix distinguishes individual mutually spatially incoherent modes, all evolving simultaneously at high speed within the beam. In the spatiospectral case, the method distinguished between mutually incoherent modes residing in the same spectral bin. Here, the method resolves between all mutually spatially incoherent fields of the beam in each temporal bin, which allows tracking their power (probability) as a function of time. We have to emphasise, however, that the method cannot directly establish the temporal coherence properties of the modes, unless a spectrometer with a resolution better than the spectral linewidth/spectral separation of the underlying spatial modes is used. The ability of the method to differentiate mutually spatially incoherent fields despite no way of separating them in time or spectrum is the added layer of flexibility enabled by applying the standard spatial tomography formalism. Whenever there is a lack of temporal or spectral resolution, the method can still spatially resolve the fields, which can be seen as a form of spontaneous spatial filtering.
We modulate the VCSEL via a pseudorandom binary sequence (PRBS) generator that is synchronised with an oscilloscope (more details of optical setup in the Supplementary Note 1). The PRBS enables repeatable temporal modulation of the VCSEL between the high (7 mA) and the low (5 mA) bias states. PRBS modulation in tandem with a synchronised oscilloscope facilitates the sequential application of several hundreds of spatial analyser states on the optical output with identical temporal modulation. We emphasize that the SLM refresh rate does not limit our temporal resolution as we collect the full temporal trace at 22 GHz sampling rategiven by the bandwidth of the photodiode usedfor each spatial analyser state displayed on the SLM.
The modulation of VCSEL disrupts the equilibrium in the cavity, which triggers dynamic changes to modes, with their relative power evolving and redistributing. Our technique allows observation of these changes for all spatial mode fields, at all delays. Figure 4 shows two temporal snapshots of the VCSEL at times t = 0.32 ns (yellow box, high bias) and t = 1.08 ns (purple box, low bias). The selected times correspond to the opposite extremes of bias applied to the VCSEL with values of 7 mA and 5 mA respectively. The overall VCSEL response (dark-grey curve, obtained from the SLM raster scan and integrated spatially for each temporal bin) closely follows the PRBS drive signal (lightgrey curve) modulating the VCSEL at 10 GHz. The remaining coloured curves in Fig. 4a are the probabilities of the spatial states that are given by the eigenvectors of the density matrices (Fig. 4h, q) in a given temporal bin.
The reconstructed spatial profiles of the first six dominant states are in Fig. 4b-g, k-p. Interestingly, during each transition from high to low bias state and vice-versa, the probability distribution of the spatial eigenstates change. The observed change is most evident on the LP21 (Fig. 4c) and LP01 (Fig. 4d) eigenstates. For high bias state, LP21 (Fig. 4c) dominates over LP01 (Fig. 4d) by about 5%. The situation is completely reversed for the low bias case, with LP01 ( Fig. 4l) having 10% higher probability compared to LP21 (Fig. 4m). This is an expected behaviour for laser diodes, with low bias favouring the low order modes, in this case the fundamental mode LP01, at the expense of higher order modes. Due to this effect, the overall spatial intensity profile of VCSEL is not constant during modulation as can be seen by comparing Fig. 4i, r obtained by SLM raster scan. The high-bias intensity profile in Fig. 4i is brighter and has more pronounced side lobes due to the higher probability of LP11 and LP21 modes compared to the low bias case in Fig. 4r which is not only dimmer but also more influenced by the LP01 mode. The observed eigenstate probability changes as a function of bias and follows similar trends to measurements of modally resolved lightcurrent curves presented in Supplementary Note 7. We also note that the SLM raster scan profiles (Fig. 4i, r) are in agreement with the probability-weighted sum of intensity profiles of all reconstructed spatial states (Fig. 4j, s), obtained using the highdimensional Stokes method. This manifests the viability of the method for temporal analysis of the beam.
Compared to the spatio-spectral analysis case, the reconstructed profiles of spatial modes appear to have distorted wavefronts. To determine the cause of this effect, we measured the variation of the VCSEL power for a fixed SLM mask (around 5%) and performed the numerical simulations of the Stokes analysis for a mixture of modes with a relative power given by the relative strength of the spectral peaks in Fig. 3a. A total of 6 modes detected in the spectral analysis (the LP31a and LP31b are not counted as they are weak) is received by the temporal detector at all times, which means that the Stokes analysis has to distinguish a mixture of 6 modes for each temporal bin (compared to 2 or 3 in spectral analysis). The numerical simulation yields similar wavefront distortion of the reconstructed modes, with increasing level of noise negatively affecting the field overlap value between the input and the reconstructed modes (see Supplementary Note 9 for detailed analysis of the noise effect on Stokes analysis). An additional source of noise can be attributed to the spatial modes of the real-world cavity being not completely orthogonal due to current injection patterns and manufacturing defects. This is evidenced by LP11 odd and even modes not having perfectly orthogonal orientation of lobes in SLM raster scans of Fig. 3g. In spectral analysis case, the effect of this imperfect orthogonality on wavefront reconstruction is not so evident because of limited number of modes in each spectral bin compared to the temporal case. We numerically simulate the effect of imperfect orthogonality on reconstructed modes in Supplementary Note 9.
The reconstructed modes and their temporally resolved probabilities/relative powers provide unique insight into the temporal dynamics within the cavity on a mode-by-mode basis. This result can be useful to explore VCSELs of varying aperture scale and shape to develop higher yield of optimum mode dynamics during fabrication. The typically used raster-scanning approaches only obtain the intensity variation of the VCSEL during modulation similar to Fig. 4i, r. While such information is useful, it does not elucidate the role of individual modes. Our tomographic approach provides not only the overall intensity profile (Fig. 4j, s), but crucially, the amplitude, phase and the relative powers of all spatial modes at each time.

Discussion
The demonstrated technique takes a major step towards enabling characterisation of an unknown beam in space along with its spectral, temporal and polarisation characteristics. The approach constitutes a simple add-on system that can analyse an arbitrary beam as it is, without access to auxiliary light sources. Moreover, as the method leverages all possible internal references, it is able to analyse even incoherent light sources where there is no suitable self-reference. Interestingly, the technique also distinguishes spatial phase, amplitude and coherence of multiple mutually incoherent arbitrary fields within a single beam, including in situations when these fields spectrally overlap or have the same time delay. By utilising the standard spatial tomography formalism, all this information is acquired in a remarkably compressed wayusing a single-pixel detector in the form of single-mode fibre and intensity-only measurements. No judicious choice of spatial basis is required as the outcome of the method is basisagnostic, which we demonstrated by using a random spatial basis throughout this manuscript. Finally, the method works for any wavelength band for which the spatial light modulator, spectrometer and photodiode exist -with most photonics labs worldwide equipped to implement the technique at almost no additional cost.
The proof-of-concept experiments explored the intricate spatiotemporal and spectral output of a modulated vertical-cavity surface-emitting laser (VCSEL) diode that represents a class of optical beams difficult to analyse using existing techniques. The lack of complete characterisation of such laser sources represents a bottleneck in a growing number of consumer, medical, automotive, sensing and imaging applications 46 . Our method unlocks information that can benefit these applications, as well as support many other areas where the full knowledge of the beam is critical, such as laser engineering 7 , optical communication in fibres and free-space 47 , quantum information 48 and remote sensing 49 . Moreover, the compatibility of the method with any intensitybased detectors opens diverse opportunities for in-depth laser diagnostics. For example, if radio-frequency spectrum analyser is used, the method can recover modally resolved relative intensity noise. The noise reduction is an important aspect of VCSEL engineering that plays a critical role in endeavours to develop Fig. 4 Spatio-temporal analysis of VCSEL (H polarisation). a Probability of the reconstructed states (coloured lines) dynamically evolves during modulation (light-grey line) and differs from the overall VCSEL response (dark-grey line). b-g and (k-p) are the first six reconstructed spatial eigenstates of density matrices (h) and (q) acquired at times t = 0.32 ns and t = 1.08 ns respectively. (i) and (r) are the total spatial intensity profiles measured by the SLM raster-scan at times t = 0.32 ns and t = 1.08 ns respectively. (j) and (s) are the total spatial intensity profiles calculated as a probability-weighted sum of intensity profiles of all reconstructed spatial states at times t = 0.32 ns and t = 1.08 ns respectively. Quantitative agreement between the reconstructed states (j, s) (tomography approach) and the raster-scan (i, r) (validation) is highlighted by the overlap integral value o evaluated in each temporal slice. The black dashed vertical lines mark the probability crossing points between eigenstate 1 and 2. For the behaviour of the system over the whole 2 ns interval see Supplementary Media 3. The V polarisation case is studied in Supplementary Note 6 and Supplementary Media 4. 100 GHz VCSELs for next-generation interconnects. Finally, based on the knowledge of the modal content, it is possible to design spatial filters that can be either used to launch one specific mode or enable observation of a target spatial mode in isolation from the rest of the system. We demonstrate these applications in Supplementary Note 10. Such spatial filter toolbox can be indispensable for establishing and removing "rogue" modes that impair the overall performance of the systemfor example, from the modulation speed or the relative intensity noise perspective.
The next-generation design of the method will aim to provide real-time feedback on the optical beam. The current configuration requires identical copies of the beam to be generated for each spatial projection measurement, which we achieve by modulating the beam with a PRBS synchronised with an oscilloscope. The total time required to display ≈ 1000 projective masks and acquire the data can be significant, currently around 30 min as detailed in Supplementary Note 9. This problem can be alleviated either by utilising multiplexed holograms 17 or by using multiplane-lightconverters 16,50,51 that both allow simultaneous interrogation of light with all spatial analyser states in parallel. The detection in such a scenario can be realised via a single-photon avalanche diode array (time) and a hyperspectral camera (spectrum). Additional improvements can be implemented in polarisation detection. Currently, our proof-of-concept system can only infer global polarisation characteristics across the whole cross-section of the beam as we do not perform state tomography on a polarisation degree of freedom. Extending the method to recover local polarisation information is possible by using the spatial and polarisation analysers in tandem, with both types of analysers implemented on respective SLMs 52 to create the required projection masks for space and polarisation degrees of freedom. Such modification would make our method suitable for the analysis of hybrid entangled states such as vector beams 53 , where the polarisation depends on the spatial position and the space/ polarisation degrees of freedom are non-separable 40 . A further serialisation of projective tomography elements (space, polarisation, energy/time) that would perform tomography of each degree of freedom is needed to extend the scope of the method to hyperentangled states in space, polarisation and time-energy 41 .

Methods
High-dimensional spatial Stokes analysis applied to each individual temporal and spectral slice. A mixture of N incoherent spatial states ψ i γ for a temporal or spectral slice denoted by γ can be mathematically represented by a density matrix ρ γ asρ where p γ i is the probability of the spatial state ψ i γ in the mixture. We emphasise that the number N and the character of the spatial states ψ i γ can be different in each spectral or temporal slice. In the following, we drop the superscript γ and elaborate on the analysis of the beam in one spectral or temporal slice. The same analysis steps are performed for all the other spectral/temporal slices. The density matrix is a positive semi-definite, Hermitian matrix with trace one. We assume that all the states in the mixture are orthogonal with respect to each other (〈ψ i |ψ j 〉 = δ ij ), which ensures the uniqueness of the density matrix. Since the density matrix is a Hermitian matrix, it can be expressed using the Gell-Mann matrices 43 where S k are equivalent to the vector elements of the Stokes vector, or, more generally, the expectation values of each Gell-Mann matrixσ k . jjσ k jj 2 F is the Frobenius norm of each Gell-Mann matrix as defined in Supplementary Note 3.
To determine the density matrix of a given mixture of states, we need to find the expectation values S k for each Gell-Mann matrix. Mathematically, this is accomplished by Experimentally, each expectation value S k is obtained via a series of projective measurements in the respective spatial basis. The details of establishing the required projections along with the experimental realisation of the projective measurements are discussed in Supplementary Note 3 and Supplementary Note 4. The reconstructed density matrix, specifically its eigenvalues and eigenvectors, encapsulates all the information about the states in the mixture and their probabilities for a given analysed spectral or temporal slice. The analysis performed over a set of spectral and temporal density matrix slices subsequently reveals the spectral and temporal evolution of multiple mutually incoherent fields in the analysed beam.
We note that the projective measurements are subject to noise, which sometimes results in the reconstructed density matrix that is not positive, semidefinite, i.e. nonphysical negative probabilities predicted. Such matrix cannot directly represent a physical mixture. We numerically find the nearest positive, semidefinite matrix by using the procedure described by Higham 45 . The effect of noise is explored and simulated in Supplementary Note 9. More details describing the working principle of the high-dimensional Stokes analysis are discussed in Supplementary Note 3.

Data availability
All data needed to evaluate the conclusions in the paper are present in the paper and Supplementary Information. Additional data related to this paper may be requested from the corresponding author.